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Abstract 

We investigate the dynamics of a Bose condensate, interacting with either re- 
pulsive or attractive forces, confined in a fully anisotropic harmonic potential. 
The (3+l)-dimensional Gross-Pitaevskii equation is integrated numerically to 
derive the collective excitation frequencies, showing a good agreement with 

those calculated with a variational technique. 
PACS numbers: 03.75.Fi, 05.30.Jp, 32.80.Pj 
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In Bose-Einstein condensation, an atomic sample confined within a magnetic trap is 
brought into a new state of matter, described by a single macroscopic quantum-mechanical 
wave function. The condensation process starts from a sample in vapour phase, and pro- 
ceeds through applications of laser cooling and evaporative cooling The condensate 
may contain between one thousand and few millions atoms. Around zero temperature or 
when the non-condensate may be neglected, the description of the macroscopic wave func- 
tion is given by the nonlinear Gross-Pitaevskii equation ||. Recently, a large attention 
within the Bose-Einstein condensation community has been concentrated on the excitation 
of condensate oscillations through a modulation of the confining potential produced by the 
magnetic trap [0-0] • The experimental and theoretical investigation of the condensate oscil- 
lation eigenmodes has provided important information on the condensate. The validity of 
the theoretical treatments has been tested, the atom coupling has been investigated in the 
different regimes, the processes responsible for the damping of the eigenmodes have been 
determined. 

In this work, we examine theoretically the condensate eigenmodes in a magnetic trap con- 
figuration leading to a very complex dynamics. Most previous analyses examined magnetic 
traps having either spherical or cylindrical symmetries, where the oscillation eigenmodes are 
determined by the trap symmetry. A different trap geometry, with different elastic constants 
along the three orthogonal axes, has been investigated in an experiment on sodium atoms, 
based on a triaxal TOP trap ||. In this letter we examine, for the first time, the collective 
frequencies associated to that magnetic trap. The calculation of the oscillation collective 
frequencies requires to solve the Gross-Pitaevskii equation (GPE) in order to determine the 
condensate spatial distribution in the ground state. Such a solution makes it possible also 
to investigate the modification of the condensate distribution following a modification of the 
magnetic trap geometry, for instance as a consequence of an oscillation of the trap confining 
potential. The GPE is here analyzed using the variational method based on the Gaussian 
ansatz for the shape of the macroscopic wave function of the condensate, applied to the 
study of the nonlinear Schrodinger equation in nonlinear optical wave propagation || and 
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for the investigation of the Bose-Einstein condensate in cylindrical Jl0[ and spherical fTT 



symmetry traps. This variational approach is here extended to the triaxial magnetic trap. 
We derive the low-energy excitation spectrum of the collective motion for both positive and 
negative scattering lengths, as a function of the atom-atom interaction strength. Our anal- 
ysis is valid for the whole range of coupling strengths, and, in the limit of large coupling, we 
recover the predictions of the hydrodynamic approach in Thomas- Fermi approximation [[| . 
We have also solved numerically the (3+l)-dimensional GPE, and derived the condensate 
mode frequencies examining the oscillation frequencies excited in the temporal evolution of 
the condensate wavefunction. The agreement between the collective mode frequencies calcu- 
lated with the approximate variational technique and those arising from the exact numerical 
integration is within a few percent. 

The macroscopic condensate wave function ip, describing an average number N of in- 
teracting bosons trapped in an external potential V at zero temperature, obeys the time- 
dependent GPE 

th^(r,t)= (^-^ + V (r) + gmf,t)\ 2 ^(r,t), (1) 

where m is the mass of the atom and the coupling constant g is given by g = Airh 2 a/m, with 
a the s-wave scattering length. The time-dependent GPE can be derived by minimizing the 
action S = J C d 3 fdt related to the Lagrangian density 

c = ihp?£ - ^fr -v^-v (f) H 2 -f H 4 . (2) 

We consider a triaxially asymmetric harmonic trapping potential of the form 

V (*0 = (A?* 2 + \W + A 2 z 2 ) , (3) 

where the adimensional constants A 2 (i = 1,2,3) are proportional to the spring constants of 
the potential along the three axes. 

Minimizing the action S within a Gaussian trial function makes it possible to obtain, for 
the condensate wave function, approximate solutions which can be derived, to a large extent, 



analytically |9|-|TT||. In the limit of weak interatomic coupling, i.e., small boson number and 



small scattering length, the choice of a Gaussian shape for the condensate is well justified 
because it describes the exact solution of the linear Schrodinger equation with harmonic 
potential. For the description of the collective dynamics of a Bose-Einstein condensate 
trapped in a spherically or cylindrically symmetric potential, it has already been shown that 
the variational technique based on Gaussian trial functions leads to reliable results even in 
the large condensate number limit [0. We take a trial function of the form 



^(f,t) = N^i — -J-— — ) 1/4 x 



(4) 

with (xi,x 2 ,x 3 ) = (x,y,z), where <7j and are time- dependent variational parameters and 
the wave function is normalized to the number N of bosons. Please note that, in order to 
describe the time evolution of the variational function, the phase factor if3i(t)x 2 is needed 

E- 

^From the Euler-Lagrange equations, the following equations of motions for the varia- 
tional parameters are derived: 

h 2 



'2 ah 2 N 



with i = 1,2,3. /,From Eq. (|5b|), we observe that the time dependence of $ is determined 
by that of <5"j. Introducing the adimensional parameters r = uot and Cj = fri/ao, where 
o-q = (h/muJo) 1 / 2 is the harmonic oscillator length, Eq. ( |5aD becomes 

—a i + X i a i = — + . 6 

dr z a? aiai<j 2 a 3 

Here, we have defined the coupling strength g = {2/ix) l / 2 Nal 'clq, proportional to the con- 
densate number N and the scattering length a, which determines the atom-atom interaction 



compared to the bare trapping potential. Equations @, with i = 1,2,3, represent a set of 
three nonlinearly coupled ordinary differential equations describing the time evolution of the 
widths of the condensate along each direction |12| : formally, they correspond to the classical 
equations of motion for a particle with coordinates Oi and total energy E = T + U, with 

T=\(al + al + at) , (7a) 
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where the dots indicate the derivative with respect to r. T plays the role of kinetic energy 
and the three terms in U, arising from the harmonic trapping potential, the kinetic pressure, 
and the atom-atom interaction, respectively, represent an effective potential energy. 

The frequencies of the low energy excitations of the condensate correspond to the small 
oscillations of the variables cr's around the equilibrium point, given by the minimum of the 
effective potential energy U in Eq. ([7b]). First, we find the equilibrium point solving the 
equations 

Ak 4 - 9— = 1, (8) 

0,jO k 

with k ranging from 1 to 3 and different from each other. 

Then, after a second-order Taylor expansion of U around the minimum, we linearize 
Eqs. (|G]) around the stationary point. The calculation of the normal mode frequencies for 
the motion of the condensate is reduced to an eigenvalue problem for the matrix A, given 
by 

d 2 U 



A 



(9) 



doidoj 

where a = (of, a®, cr®) stands for the solution of Eq. @. In the Thomas-Fermi limit, i.e., for 
g ^> 1, the right-hand side in Eq. (||), related to the kinetic pressure, can be neglected and an 
analytic expression for the equilibrium point is obtained. In this case, the eigenfrequencies 
uj for the collective motion, in units of ujq, are given by the solutions of the equation 



= u 6 - 3 (A? + X\ + X 2 3 

+8 (Xl\l + XjXl + XjX?) oo 2 - 20XlX 2 2 Xl (10) 

The same equation, obtained within a hydro dynamic approach, is reported in Ref. H. In 
the opposite limit, for g = 0, we recover the well-known result Ui = 2Aj, with i — 1,2, 3, in 
units of u , for the normal modes in the non- interacting case. 

The numerical integration of the time dependent GPE, in Eq. (|l|), was done using a 
modified split operator technique, adapted to the integration of a Schrodinger equation. We 
wrote Eq. (jl]) in the form 

ih^{r,t) = (H x (f,t) + H y (f,t) + H z (f,t))^(f,t), (11) 



where 



h 2 d 2 .. l 



Hi(r,t) = -^^f + Vte) + -gmm 2 . (12) 

The idea is to split the full Hamiltonian in three sub-Hamiltonians, so that at each time we 
have to write the Laplacian with respect to one coordinate only, leading to the solution of a 
tridiagonal system, and to huge savings in computer memory. The splitting is carried out so 
that the commutators are exact up to the order S 2 included. Equation (|ll]) was integrated 
using the scheme (5 is the integration time step, and Ai(t) = i5Hi(r,t)/h) 

' ;i-A z (t)/2)x 



1 + A,(t)/2 

T -±_ { l-A,(t ) /2H.(f^ (13) 

There is obviously a problem with the nonlinear term g\ip(r, t) | 2 , because we should really use 
a ip somehow averaged over the time step 5, not a ip evaluated at the beginning of the time 
step. To circumvent this problem, we used a sort of predictor corrector step. Each integration 
step is really done in two times: going from the time t to the time t + 8, the first time we 
used ip(r,t) in the nonlinear term, obtaining a "predicted" ip(r,t + 5); we then repeated 
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the integration step, starting again from ip(f,t), but using | \ ip(r,t) + ip(r,t + 5) J in the 
nonlinear term. In Fig. [l| we show a typical time evolution obtained from the integration 
of the GPE. We plotted here N' 1 J \ip(r, t)\ 2 dy dz as funct ion of t and x. The initial wave 
function we started from has the form of Eq. (ffl), with <7; = 1 and $ = for i — 1, 2, 3, i.e., 
fairly far from equilibrium. 

The Gaussian wave function is not an exact steady-state solution for the GPE, as 
also stated in [13[]. The small deviations from the actual stationary solution excite small- 
amplitude oscillations in the condensate distribution. We have spectrally analyzed the time 
dependence of the condensate widths along the three spatial directions, in order to derive 
the oscillation frequencies. This analysis has been done for different values of g, for both 
positive and negative scattering lengths. A typical spectral density is shown in Fig. § (top). 

The general solution of Eq. (§) for the stationary state and of Eq. the eigenvalue 
problem for the collective mode frequencies, was obtained numerically. This solution is shown 
in Fig. H] (bottom), where the scaled widths <7j of the stationary condensate wave-function, 
solid lines, and the normal mode frequencies cjj, dashed lines, are shown as functions of 
the coupling strength g in case of a positive scattering length a. Figure ^ (top) shows 
the same quantities for a negative scattering length a. On the top axes of the graphs in 
Figs. @ (bottom) and § (top), the number N of atoms in the condensate is reported, for 
scattering lengths and trapping potential parameters fixed to the values listed below. The 
markers in Figs. [| and ^| represent the mode frequencies calculated from the numerical 
solution of the GPE. The results in Fig. |2] are obtained for parameters corresponding to 
the triaxially anisotropic trap used in the experiment in Ref. ||, performed at NIST on Na 
atoms, e.g., tu = 2ix x 234 s _1 , Ai = 1, A 2 = V2 and A 3 = 2. The Na scattering length 
a = 2.75 nm is used in the calculation. The equilibrium solutions for the condensate widths 
(jj are increasing functions of the condensate number N and, asymptotically, scale as iV 1 / 5 . 
The mode frequencies decrease with N monotonically from the non-interacting values to the 
asymptotic values, indicated in the figure by the solid horizontal segments. 
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Figure ^| (top) shows the solution for a negative scattering length a = —1.4 nm, appro- 
priate for 7 Li atoms. The parameters characterizing the trapping potential are the same as 
in Fig. |2|. For small condensate numbers we recover again the non-interacting limit. With 
increasing N, the atom-atom interaction strongly modifies the dynamics of the condensate, 
in a different way with respect to the case of positive scattering length. Notice the differ- 
ent horizontal scales in Figs. ^| and It is worth noting that the variational formalism, 
being applicable for any condensate number N, makes it possible to investigate the case of 
negative scattering lengths. In this case, the large N limit of the Thomas-Fermi approx- 
imation and the hydrodynamic approach do not apply, because the number of condensed 
atoms cannot exceed a critical value [§]]. For negative coupling constant g, Eqs. (||) present 
either two solutions for the stationary width <jj of the condensate along each direction, or 
no solution, depending on the absolute value of g. This is shown in Fig. |3|, where, for each 
<jj, a two-branched solution is displayed (solid lines). The upper and lower branches of the 
o"j curves correspond to the stable and the unstable solution, respectively. As functions of 
\g\, the stable stationary widths shrink, while the unstable ones, originating from zero as a 
limit for \g\ — > 0, increase. The two branches of each cij solution merge at a critical value of 
the coupling strength. Beyond this value, no solution to Eqs. (]8|) exists. The dashed lines in 
Fig. § show the frequencies of the collective motion around the stable solution. At the criti- 
cal point, the two highest frequencies diverge and the lowest one falls to zero. Similar results 
have been obtained in Ref. jl0|| for cylindrically symmetric condensates. This theoretical 



results agree well with the numerical integration of the GPE. 

However, in the case of negative scattering lengths the collapse of the condensate implies 
that the GPE can be numerically integrated only until the collapse takes place. As the 
number of particles in the condensate increases and approaches the critical value, the lifetime 
of the condensate decreases and, eventually, becomes so short that we can no longer observe 
any oscillatory dynamics. For the Li parameters, the collapse lifetime is shown in Fig. |3] as 
function of N, for iV approaching the critical value. 

In conclusion, we have performed an analysis of the dynamics of a Bose-Einstein con- 
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densate trapped in a triaxially asymmetric potential. Condensation was achieved at NIST 
in a fully asymmetric TOP trap ||, while a specific study of the condensate dynamics in 
such a geometry was lacking. At present, no experimental data concerning the collective 
motion of the condensate are available from the NIST experiment. Our predictions might be 
confirmed by forthcoming experimental results. Within the different point of view started 



by Dalfovo et al. ||14|| , the study of the nonlinear coupling between the condensate modes 
provides a stronger test of the condensate characteristics. Moreover, nonlinear processes 
could be useful to perform nonlinear spectroscopy or nonlinear optics in the condensate. As 
well known in the framework of nonlinear dynamics, more complex scenarios are produced 
when the number of degrees of freedom is increased, as in the (3+l)-dimensional evolution 
in a triaxial magnetic trap. 

This project was supported by the Italian INFM through a PRA project on Bose-Einstein 
Condensation and by the CNR through a Progetto Integrato. The authors wish to thank 
L. Reatto for a careful reading of the manuscript. 
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FIGURES 

FIG. 1. A typical result of the numerical integration of the GPE: plot of N^ 1 J \ip(r, t)\ 2 dy dz 
vs t and x. We used the Na parameters for the integration (see text), and g = 3.335. 

FIG. 2. Top: Typical spectral densities from the numerical integration of the GPE: we show 
the Fourier transform of the time evolution of the condensate widths along the three axes, for the 
Na case and g = 0.664. The arrows show the location of the eigenmode frequencies expected from 
the variational analysis. Bottom: Results for the variational calculation of the scaled oscillation 
frequencies LOi, in units of loq (dashed lines), and scaled widths <7j, in units of ao ( continuous lines), 
as a function of g (lower scale) and N (upper scale). The horizontal segments on the right represent 
the hydrodynamic limits. The markers indicate the oscillation frequencies derived from numerical 
integration of the GPE. Triaxial TOP parameters as in (|| and scattering length a = 2.75 nm. 

FIG. 3. Top: As in Fig. [2|, except scattering length a = —1.4 nm. Bottom: time before we 
observe the condensate collapse, as a function of N, for the same parameters of the top figure. The 
initial wave function has the form of Eq. (|j), with the Oi equal to the variational ones, and the 
A = 0. 
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